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Abstract 

A recently developed method [A. Shelankov and M. Ozana, Phys. Rev. B 61, 
7077 (2000)] is applied to investigate d-wave superconductors in the vicinity 
of (rough) surfaces. While this method allows the incorporation of arbitrary 
interfaces into the quasiclassical technique, we discuss, as examples, diffusive 
surfaces and boundaries with small tilted mirrors (facets). The properties 
of the surface enter via the scattering matrix in the boundary condition for 
the quasiclassical Green's function. The diffusive surface is described by an 
ensemble of random scattering matrices. We find that the fluctuations of 
the density of states around the average are small; the zero bias conductance 
peak broadens with increasing disorder. The faceted surface is described in 
the model where the scattering matrix couples m in- and m out-trajectories 
(m > 2). No zero bias conductance peak is found for [100] surfaces; the 
relation to the model of Fogelstrom et al. [Phys. Rev. Lett. 79, 281 (1997)] 
is discussed. 
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I. INTRODUCTION 



The (ia;2-y2 symmetry of the order parameter (OP) in high temperature (high-T(.) super- 
conductors is nowadays well established by various phase-sensitive experiments, the most 
striking being the observation of a half flux quantum in a tricrystal geometry by Tsui 
et al. [||. Also the zero-bias conductance peak (ZBCP) found in different tunnel exper- 
iments on [110]-oriented boundaries of YBa2Cu307_5 gave a strong hint for d-wave 
symmetry as was first pointed out by Hu 0|. However, in more recent experiments, several 
puzzling questions arose. For example a splitting of the ZBCP was observed and ZBCP's 
were reported even for [100]-oriented junctions Disorder effects were examined as well 



in experiment where the disorder of the junctions was increased by ion irradiation |jTOl- A 
decreasing height of the ZBCP was observed whereas the width stayed constant. 

The most successful approach to treat such inhomogeneous problems is the theory of 
quasiclassical Green's functions |TT|-|T3|. The properties of surfaces or interfaces are included 
by effective boundary condition. For a specular surface the Green's function has to be 
continuous on a classical trajectory (see Fig. |l|). In this most simple model the pair-breaking 
effect of surfaces as well as the existence of a ZBCP can be explained: If quasiparticles are 
scattered to branches with a different OP (a 7^ ) the pairing is suppressed and bound 
states can occur due to Andreev scattering. If a quasiparticle is scattered with a sign change 
of the OP the ZBCP exists, too. In this framework it was also shown 0,11^] that the 



splitting of the ZBCP is in agreement with the existence of a subdominant order parameter 
{dx^-y2 +idxy/ s) in the vicinity of the surface. The generalization of the boundary conditions 
to a specular interface was carried out by Zaitsev ||16| . In this situation the Green's function 



on four trajectories must match at the interface, which leads to quite complicated non-linear 
boundary conditions. 

As surface roughness is present in experiments, models were developed to include disorder 
in the theory. One possibility, first suggested by Ovchinnikov |T^, is to use the boundary 
conditions for the specular situation and cover the surface with a thin dirty layer where 
equations for the dirty case must be applied. In numerical studies of boundary problems 
similar techniques were used |T^,[T^]. Also a scattering matrix approach, which is related to 
the randomly rippled wall model [^, was applied to disordered surfaces [Q; the solution 



was given in a Born-like approximation. Except for models with unitary scatterers [^, all 
calculations show a broadening of the ZBCP due to disorder. 

Often the surface roughness is only present on scales much smaller than the coherence 
length. In this case, the boundary conditions for the quasiclassical Green's function can be 
formulated in terms of the scattering matrix (^-matrix), as it has been recently suggested 
in A rough interface does not conserve the momentum parallel to the surface, and, 

therefore, it couples waves {i.e. quasiclassical trajectories) with different propagation direc- 
tion. The microscopic structure of the interface enters the theory via the S'-matrix. In the 
absence of detailed knowledge about the microscopic structure of the surface, the S'-matrix 
has to be taken phenomenological input. 

Due to the very short in-plane coherence length in high-Tc superconductors (^0 ~ 20A) 
structures on larger scales occur as well. Then the translational invariance parallel to the 
surface is lost on this length scale and a full 2D treatment of the problem is necessary, in 
which the trajectories are considered individually (see Fig. |^). For example, in cuprates 
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facets with typical dimensions of 10 — lOOnm are present at interfaces p3[. This leads to 



the existence of a ZBCP even for [100] tunnel junctions as was pointed out by Fogelstrom 
et al. 



In our study we use the scattering matrix approach presented in |2J] to describe surfaces 
with microscopic roughness. We will discuss two kinds of surfaces: First, we study a micro- 
scopically disordered surface, which is described by random matrices. In contrast to earlier 
calculations we are able to consider individual realizations of the disorder; we examine aver- 
aged quantities as well as fluctuations. Afterwards we focus on an S'-matrix that describes 
a surface with tilted mirrors, where few trajectories are connected at the surface. 

In the following section we will briefly introduce the quasiclassical theory for supercon- 
ductors. Subsequently we will present the boundary conditions in the form most suitable 
for our purpose and discuss the general properties. In section |IT1| we derive the S-matrices 
for different situations and present the related results. We conclude with a discussion of our 
results and compare them to other approaches. 

II. METHOD 

A. Theory of quasiclassical Green's functions 

In our studies of boundary effects of unconventional superconductors we use the theory 
of quasiclassical Green's functions introduced by Eilenberger |]TT|. This approximation of 
Gorkov's theory is valid in the quasiclassical limit {271 /k-p <^ 0- Several review articles 



have been written on this subject, e.g. by Schmid |12[ or Larkin and Ovchinnikov [|13l. The 
quasiclassical propagator in Nambu space 

^(^,R,kp)= (jt /t) (1) 
is determined by the Eilenberger equation 

(u + -WY ■ A(R)) fa - A(R, kp), ^(cu, R, kp)] + ^^vf • V)^(cu, R, kp) = (2) 
where fj represent the Pauli matrices and 

A(R,kp)=(_^, fj. (3) 

Additionally a normalization condition is needed to obtain the physical solution of the 
equation 

^2(cu,R,kp) = i. (4) 
The OP must obey the self-consistency equation 

— (\/(kp,k^)^^(c.,R,k^)>^, . (5) 
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Here (. . . indicates the average over the Fermi surface. In thermal equihbrium the 
Keldysh propagator is given by the the advanced and retarded propagator g^/^ 

g^ = (g^ - g^) tanh(cu/2A;BT) with g^^^ ■.tu^uj±iO+. (6) 

For simphcity we make some further assumptions concerning the microscopic properties: For 
the interaction we choose l^(kF, kp) = V cos[2{(p — a)] cos[2{ip' — a)] which generates a 
(i-wave OP with orientation a (see Fig. 

A(R,kF) = A(R)cos[2(v9-a)]. (7) 

In addition we assume an isotropic two-dimensional model with a spherical Fermi surface. 
After the determination of the OP all physical properties can be calculated from the qua- 
siclassical Green's function. For example, the angle-resolved local density of states (DOS) 
reads 

Xiu, R, kp) = AAoRe [g'^iu, R, kp)] , (8) 

where Ao is the normal state DOS. In many cases the knowledge of the angle-averaged DOS 
is sufficient 

^^{uJ,R) = {^^{cu,RM)K■ (9) 

The DOS at the surface can directly be measured via the differential conductance G = dl/dV 
for normal-metal-insulator-superconductor tunnel-junctions. For T — it is given by [jl8| 



G{V) = e^A{vF,.TikF)^^{eV,x = OM))ik,^.>o), (10) 
where the transmission probability is chosen as 

T(^) =t2sin2((^) < 1, (11) 

and A is the area of the contact. The current- density can be calculated from the Keldysh 
Green's function via 

j(R) = -eATo J ^ (vpTr[f3^^(^, R, kp)]>^^ . (12) 

It has been shown that the decomposition introduced by Maki and Schopohl |2^ is suitable 
for the numerical integration of the Eilenberger equation as well as for analytical consider- 
ations (see [II B|) 



1 + ab —2a 



ab V 26 -{l + ab)l ' 



(13) 



Considering the physical meaning, the functions a{uj, R, kp) and b{uj, R, kp) are closely re- 
lated to the particle and hole amplitudes in the Andreev equation as was discussed in detail 
23. With this construction the normalization condition is obeyed automatically. By 



m 
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putting in this decomposition in Eq. (0) it can be seen that the functions a and h are given 
by the equations 



ih{wY ■ V)a = A*a^ - 2uja + A, (14) 
-i^(vF ■ V)b = A6^ - 2ujb + A*. (15) 

These equations can be solved on classical trajectories labeled by the Fermi wave vector kp. 
For each direction vp one has to integrate two ordinary differential equations in order to 
construct the full propagator. 

The Matsubara technique can be used as well to calculate the OP 

AiR,kp) = -kj,Tm J2 (l"(kF,k^)/^(u;„,R,k^)) , (16) 

\LlI„ I <LUc 

and the current- density 

j(R) = -eAfokBTm (vFTr[f3^^(a;„, R, kF)])^^^ . (17) 

The energy integrals turn to sums over discrete Matsubara frequencies Un = k^TT:{2n + 1) 
and the Matsubara propagator is determined by the relation 

^^/^(c.,R,kF) = ^^(^„,R,kF)|^^^_^„^. (18) 

One crucial point for investigating the effects of boundaries is still missing. As the qua- 
siclassical condition does not apply in the vicinity of surfaces and interfaces we have to 
treat the scattering of quasiparticles by effective boundary conditions. The properties of the 
boundary enter the calculations only at this point. 



B. Boundary Conditions 

In our work we use the general theory recently derived in ||2^. The starting point is 
the Andreev-like equation for the particle- and hole-like amplitudes which factorize the 
Eilenberger Green's function in Eq. (see for details). In this approach, it is possible 



to consider roughness that occurs on length scales much smaller than the coherence length. 
All information on the microscopic shape of the boundary is provided by the scattering 
amplitudes from the in-trajectories (/cpa; < 0) to the out-trajectories {k^l > 0); they are 
gathered in the scattering matrix S. 

For simplicity we consider only a finite number n of discrete in- and out-trajectories 



j^m/out _^ j^m/out^ ^ _ 2 . . . with equidistant angles. Following the boundary condi- 
tions are determined using the functions 



Ai{(3) = det[l - SaS^Sf], (19) 
Bi{a) = det[l - Sa'j'S^b] (20) 
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with the diagonal n x n-matrices 

d = diag{ai, . . . , a„}, a" = diag{ai, . . . , ai_i, a, a^+i, . . . }, 

6 = diag{6i, . . . ,6„}, 6f = diag{6i, . . . , /3, . . . }, (21) 

and flj = a{uj, x = 0, kp J, 6j = b{u!, x = 0, kp*^*). 

The solutions of y4/(/9) = and Bi{a) = provide the boundary conditions 

A{(3o)=0 a(a;,x = 0,k°';f) = l, (22) 

Po 

Bi{ao)=0 6(o;,a; = 0,ki?;) = — . (23) 

As the determinant is a linear function of each of the matrix elements the functions Ai{f]) 
and Bi{a) are linear in /? and a. We can solve the boundary condition by calculating Ai{P) 
and Bi{a) for two arbitrary values of /3 and a; for (3 = 0,1 and a = 0, 1 we obtain 

a(^,x = 0,k-*) = l-^, (24) 

6(^,a; = 0,k|?,) = l-||||. (25) 

With the boundary condition the Green's function can be calculated: At first the integration 
of Eq. (|14|) on the in- and of Eq. (|1^) on the out-trajectories is performed starting from the 



known bulk values \24 



a{u, X oo, ki?) = ^^=M_=, (26) 



,out\ 



A* (\(° , 

h{uj, X ^ oo, kD = oo^ ^ ! _ (27) 



towards the boundary (Aqq: bulk OP). Then the boundary conditions must be applied to 
get the a's on the out- and the 6's on the in-trajectories at the boundary and the succeeding 
integration on these trajectories provides the missing a's and fe's. 

The properties of the boundaries enter only via the 5- matrix. The value \Sij\'^ is the 
probability of scattering from kp^- to kp*^*. We choose the numbering of the trajectories so 
that S = 1 reproduces the specular case. Due to current conservation S must be unitary 

SS^ = 1. (28) 

With a suitable choice of S arbitrary physical situations can be treated by this technique. 



Some examples are presented in chapter |T] 



We are also able to connect basic symmetries of the physical situation with transfor- 
mation properties of the scattering matrix. The symmetry operations for the mirror and 
the time-reversal symmetry are illustrated in Fig. |^. The mirror symmetry of the surface 
[y —y) is given by the transformation 

S' = TST, with = 5(^n+l-^),j. (29) 

The time-reversal symmetry operation is represented by 

S' = TS^T. (30) 
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III. RESULTS FOR DIFFERENT SURFACES 



Since we are discussing different types of rougliness which occur in experiments we have 
to find adequate scattering matrices for each situation. As the unitarity condition (^) must 
be obeyed we represent S by the relation 

S = exp{iH} with H = (31) 



In the subsequent sections [111 A| and [III Ij| we present S-matrices for random surfaces as 



well as for surfaces with small tilted mirrors and physical properties such as the OP and 
the DOS in the vicinity of a surface. The retarded Green 's-f unction should be evaluated 
at u uj + i6] for numerical purposes we keep 6 finite. For the calculation of the DOS we 
choose 6 = 0.02kBTc and n = 200; we checked that the results do not change for larger n. 



A. Random surface 

In this section we search for a scattering matrix that can describe random surfaces. To 
take into account the statistical properties of the surfaces we choose a random matrix in this 
approach; therefore the Hermitian matrix H is assumed to be a member of the Gaussian 



unitary ensemble (GUE) as was also suggested by Yamada et al. |^ 

{H^,) = 0, 



r 



(32) 



The brackets (...) denote the ensemble average of the disorder. The roughness of the surface 
can be varied by the parameter r. The factor 1/n in the correlator of H ensures that the 
whole procedure does not depend on the number of channels for n — oo. 

The averaged scattering probabilities {\Sij\'^) have a simple behavior. For r = only 
1 5*42 p = 1 are finite and all other elements are zero. If r is increased the diagonal elements 
(responsible for specular reflection) are reduced to {\Sii\'^) = |m(t)P < 1 and the off-diagonal 
elements become finite {\Si-^j\^) = \v{t)\'^ < 1/n (see Fig. ^) with 

\u{r)\' + {n - 1Mt)\' = 1. (33) 

Therefore, the averaged properties of S are fully determined by the probability for specular 
reflection |mP and we can use it as a measure for the disorder of the surface; its relation to 
the parameter r is shown in Fig. ^. For r small enough, the reflection is partially specular. 
When r > 2, the scattering becomes isotropic since \u\^ ~ |f p ~ 1/n. We call this situation 
the diffusive limit. 

We apply a random matrix S to calculate the OP and the DOS in the vicinity of a 
disordered surface. We study several (up to 50) realizations of the S'-matrix individually. 
We find that the fluctuations are small as can be seen in Fig. Therefore, it is meaningful 
to consider the averaged quantities (A) and (G). The results for a = 0°,45° and different 
roughness values are shown in figs. |^ and |] for the OP and for the differential conductance 
in figs. ^ and |TD|, where the normal state resistance is used 

RJ,' = e^AMot^VF^. (34) 

OTT 
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Additionally the angle-dependent DOS is presented in Fig. [TT| for the medium roughness 
r = 0.4. 

We point out some interesting features of the model: Disorder leads to a suppression and 
broadening of the ZBCP for a = 45°. By comparison with experimental data [l^ realistic 
results can be achieved by 0.8 ^ t < 2. In our model with disorder, for a = 0° no ZBCP 
occurs. In the angle-resolved DOS no splitting of any bound states due to disorder is seen. 
In the diffusive limit the OP reaches an almost universal curve independent of the surface 
orientation a; the conductance becomes flat and is of the order of the normal state value for 
all energies and surface orientations. 



B. Surface with small tilted mirrors 

With this method also other roughness types can be studied. In this section we consider 
surfaces which only connect few trajectories. For example one can think of small (compared 
to the coherence length) mirrors with distinct orientations. Each of the mirrors contributes to 
reflection so that the surface acts as a beam-splitter. For the mirror orientations 6 = 180°//m 
with integer / = — m + 1, . . . , m — 1, this can be achieved by the choice 

m—l 

H = rH with Hij = ^ 5\i-j\,in/m- (35) 

1=1 

In the further calculation we choose n to be a multiple of m; this leads to a simpler form of 
S but has no physical implication for n ^ m. In this case H considered as a block m x m 
matrix has zeros on the diagonal, and all other elements are equal to l„/m {{n/m) x (n/m) 
unity matrix). Since 



the exponential representation for the S-matrix in Eq. ( ^T]) becomes 



where 



= (m - 1)1 + (m - 2)H 



S = M(r)l + t;(r)H, 



\u\^ + (m - l)\v\'^ = 1, 

uv* + u*v+ (m-2)|t;p = 



(36) 

(37) 

(38) 
(39) 



as required by the unitarity of the S"- matrix. For fixed r the value of |mP is the probability 
for specular reflection, whereas |f p is the probability for the scattering on one of the tilted 
mirrors. The scattering matrix can be calculated explicitly and the amplitudes are given by 

M(r) = ^(m-l + e^'^n, (40) 
m 

v{t) = ^(e^"^" - 1). (41) 
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Therefore the probabihty for specular reflection (mirror with 6 = 0) is given by 

\u{t)\^ = -^[{m - 1)2 + 1 + 2(m - 1) cos(mr)]. (42) 

In this model, each incoming trajectory is split into m outgoing trajectories. Altogether, 
there are m in-trajectories which are scattered into the same out-states. In the terminology 



of [|^, this corresponds to a knot with m in- and m out-trajectories (see Fig. |T^). 

We study the simplest case m = 2 with three mirrors with orientation 6 = 0°,±45°, 
where two in- and two out-trajectories are coupled. The S'-matrix has the form 

with u = uln/2 and v = vln/2', the functions u and v are given by 

u{t) = cos(r), v{t) = 'isin(r). (44) 

For a = 0° and Iwp < 1 finite energy bound states occur; their energies move to zero with 
decreasing weight of specular reflection. If \u\'^ = (r = 7r/2), i.e. no specular reflection 
is present, the bound states reach zero energy. For a = 45° this is reversed: non-specular 
reflection leads to a splitting of the zero energy bound states, which grows with decreasing 
\u\'^ (Fig. |13[), as was discussed qualitatively in [Q. This model is in some aspects similar 
to a Josephson-contact with a specular interface, where also a splitting of the zero energy 
bound state can be observed. 

The case m > 2, however, is more complex and cannot be mapped to any studied 
model. We consider the case m = 3. The five mirrors have the relative orientations 6 = 
0°, ±30°, ±60° and three in- and three out-trajectories are connected; the S'-matrix is given 
by 

/ U V V \ 

(45) 

with u = uln/3 and v = vln/3- As can be seen in Eq. (|4^ ) for m = 3 it is only possible to 
choose the specular scattering probability in the interval |mP G [1/9, 1]. For a = 0° non- 
specular scattering leads to finite energy bound states which are moving to lower energies 
with decreasing |up; but here zero energy is not reached. In the case of a = 45° for \u\'^ < 1 
one part of the zero energy bound state splits to finite energies, whereas another part stays 




at zero energy with reduced spectral weight (Fig. |T4D. In figs. |T5| we also present the case 
m = 4 with 7 different mirrors for the surface orientations a = 0°,45°. 

Summing up our observations for the a = 45° case small mirrors reduce the spectral 
weight of the zero energy bound states as a part splits to finite energies. For the surface 
orientation « = 0° in general (except some particular situations) no zero energy bound state 
is produced. 

Finally we point out that a larger class of S'-matrices can be used to describe such surfaces 

m—l 

Hij = ^ n5\i^j\^ln/m- (46) 
1=1 
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This model has almost the same properties as (|35D, however, the scattering probabilities for 
each mirror differ. Moreover it is possible to combine the these mirrors with disorder just 
by multiplying the related S'-matrices and averaging as in section [111 A| . 



IV. DISCUSSION AND CONCLUSION 

In the present paper, the scattering matrix approach has been applied to describe d-wave 
superconductors in the vicinity of rough surfaces. Two physical situations are examined: (i) 
a surface with partially diffusive reflection described by random scattering matrices; (ii) a 
surface with small tilted mirrors (facets) where the reflected wave is a coherent mixture of 
waves propagating in several directions. 

First, for the diffusive surface it appears that our results are very similar to those found 
in other approaches such as the randomly rippled wall model or the thin dirty layer |]18 



In contrast to those calculations we treated the disorder by direct sampling. We find that 
the deviation of the DOS from the average is rather small. Therefore, our calculations con- 
firm the validity of averaging procedures used in the afore-mentioned papers. In particular 
the broadening of the ZBCP due to increasing disorder is no artifact of the approximate 
averaging. (Other models exist that show no such broadening pD| , p7|] .) 

Second, we studied a surface with tilted mirrors. In contrast to the model for large 



facets examined in , our model describes a surface where faceting occurs on a scale small 
compared to the coherence length. These two models provide qualitatively different results: 
our model in general shows no ZBCP for a = 0, whereas large facets lead to a ZBCP for 
each surface orientation. 

In experiments on high-Tc materials, it cannot be excluded that roughness on a scale of 
the coherence length or larger occurs, which is beyond the model used in the current paper. 



This might be the case in the experiment described in [10|, where the width of the ZBCP 



is constant with varying disorder. In some experiments a ZBCP for a = is observed, 



too, consistent with the model for large facets |15 . 
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FIG. 1. Surface effect in real space (left) an k-space (right). 
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FIG. 2. Effect of roughness on small (left) and large (right). 
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FIG. 3. Two trajectories (solid and dashed line) arc plotted which are connected via symmetry 
transformation for mirror (left) and time-reversal symmetry (right). 
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FIG. 4. Mean scattering probability from the in-trajectories [i + I) to a fixed out-trajectory 
z = n/2 for different values of r. The specular contribution is reduced by increasing r. 
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FIG. 5. Specular scattering weight as a function of r. The diffusive limit is reached for r > 2. 
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FIG. 6. Conductance for a = 45° and r = 0.8, 2; the typical fluctuations are confined by the 
dashed lines. 
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FIG. 7. Averaged order parameter for a = 0° and r = 0, 0.08, 0.4, 0.8, 2, 4. 
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FIG. 8. Averaged order parameter for a = 45° and r = 0, 0.08, 0.4, 0.8, 2, 4. 
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FIG. 9. Averaged differential conductance for a = 0° and r = 0, 0.08, 0.4, 0.8, 2, 4. 



15 



2 



~\ a = 45° T = 



T = 0.08 




0.2 0.4 0.6 0.8 1 1.2 

FIG. 10. Averaged differential conductance for a = 45° and r = 0, 0.08, 0.4, 0.8, 2, 4. 




0.2 0.4 0.6 0.8 1 1.2 

FIG. 11. Averaged angle-resolved DOS for a = 45° and r = 0.4. 




FIG. 12. For m = 3 each of the three in-trajectories (solid lines) contributes to the same three 
out-trajectories (dashed lines). 
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FIG. 14. Differential conductance for m = 3 and a = 0°,45°; the specular scattering weight is 
varied. 
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FIG. 15. Differential conductance for m = 4 and a = 0°,45°; the specular scattering weight is 
varied. 
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